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Abstract 

This paper develops a first-order system least-squares (FOSLS) formulation for equations 
of two-phase flow. The main goal is to show that this discretization, along with numerical 
techniques such as nested iteration, algebraic multigrid, and adaptive local refinement, 
can be used to solve these types of complex fluid flow problems. In addition, from an 
energetic variational approach, it can be shown that an important quantity to preserve 
in a given simulation is the energy law. We discuss the energy law and inherent structure 
for two-phase flow using the Allen-Cahn interface model and indicate how it is related 
to other complex fluid models, such as magnetohydrodynamics. Finally, we show that, 
using the FOSLS framework, one can still satisfy the appropriate energy law globally 
while using well-known numerical techniques. 

Keywords: multiphase flow, energetic variational approach, algebraic multigrid, 
first-order system least squares, nested iteration 



1. Introduction 

Complex fluids involve multi-physics and multi-scale phenomena that require ad- 
vanced techniques in scientific computing in order to be solved efficiently. The goal of 
this paper is to present numerical techniques that have been well tested and shown to 
give accurate solutions with a reasonable amount of computational cost, while also pre- 
serving the global energy of the system. While the systems can be described with partial 
differential equations (PDEs), in this case a time-dependent and nonlinear system of 
equations, an underlying energy law also describes the physics. It is from these laws that 
the PDE is derived. Therefore, approximate solutions to complex fluid problems should 
adequately approximate this energy law. 

Here, we look at one type of complex fluid problem, two-phase flow. The conventional 
sharp interface description of the mixture of two Newtonian fluids with density, pi, and 
viscosity, /i,-, assumes the mixture occupies the overall domain, Q = fii U SI2, as in figure 
[T] The interface between two materials is the free moving interface, T t . 
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Figure 1: Fluid 1 is in region f2i and fluid 2 is region S72. Tt represents the mixing interface. 



The overall system includes the Navier-Stokes equation in each f^, where i = 1 or 2: 



V dt 



Vu 4 = 0. 



Here, u^, pi, pi, and /Ltj are the velocities, pressures, densities, and viscosities of each 
fluid, respectively. Two types of boundary conditions are possible. Kinematic boundary 
conditions yield Uj-n = V n , where n is the "directional" normal of T t and V n is the normal 
speed of I\. Traction free boundary conditions (also known as force balance conditions) 
satisfy [rj] • n = surface force on T t , where [r^ is the jump of the stress tensor, 

n = + pil, £. = - . 

Remark Assume u 2 |ar2 = 0. With the traction free boundary condition, this implies 

r ( nSO = 0. 

In case the surface force is attributed to the surface tension, i.e., F gur f ace = kHn, where 
H is the mean curvature of T t and k is the surface tension parameter, then the whole 
system possesses the underlying energy law, 

a {* r i ' 2 



dt 



ip 4 |u 4 | 2 dx + fcArea(r t )j - - \Y,f £ Mi|Vu;| 2 dxj . 

In this paper, we take a different approach, called the diffusive interface method 
(phase- field method) . The label function, </>(x, t) , is introduced such that 



0(x,i) 



T in fluid 1, 
•1 in fluid 2. 
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The interaction between the two materials is reflected through the "mixing" energy, or 
the Ginzburg-Landau type of energy: 

f W(cb)dx = [ I|V^| 2 + -^(0 2 -l) 2 dx, 

J- O J- f2 

where J Q ^|V0| 2 <ix represents the tendency of the mixture to be homogeneous. In other 

words, it is the "phillic" interaction. In addition, ^ (</> 2 — l) 2 cfec represents the 
"phobic" interaction, or the tendency of the mixture to be separated. The parameter e 
reflects the competition between the two opposite tendencies. 

Due to the nature of this empirical energy, it is present in many theories in physics, 
such as superconductivity Q3 [2] and phase transitions [3J 0j. It can be shown that 
<f> — > ±1 almost everywhere as e — > 0, thus representing the separation of the two im- 
miscible materials. Moreover, the thermodynamics and patterns of the interface, T t , 
can be captured by the level sets of <f>, for instance, {x|<^(x, t) = 0}. In addition, the 
energy e J Q W((f>) — > coArea(T t ), which is equivalent to the surface tension case in the 
sharp interface formulation 5 . It is worth pointing out that, while one can view the 
diffusive interface formula as being the regularizing approximation of the sharp interface 
formulation, it is more physical to view the sharp interface formulation as being the ideal 
approximation of the diffusive interface method. Therefore, it can be understood that 
the phase field description really captures the microscopic interaction of the mixtures. 

In the examples presented here, we consider problems where we have one or more 
bubbles of one phase interacting in a region of a fluid of the other phase. From this 
system, the continuous energy laws are derived as is the continuous system of PDEs. We 
then show that these laws are adequately preserved numerically. 

To solve the system, we use the first-order system least-squares (FOSLS) finite ele- 
ment approach. This method has been well-tested on many different types of PDEs and, 
specifically, on several complex fluid systems [H [71 |H1 El HH] > including those of magneto- 
hydrodynamics (MHD) [H1[T2|. We show that use of the FOSLS method allows both an 
accurate solution of the system and preservation of the energy laws. Moreover, because 
we are using this method, we have the advantage of being able to use other numerical 
techniques, such as nested iteration and adaptive local refinement, which greatly improve 
the efficiency of solving these problems. 

This paper is structured as follows. In Section[2j we describe the two-phase model 
and present the equations to be solved. In section |3j we discuss the FOSLS formulation 
that is used and show how it presents a well-posed discrete problem. Section [4] introduces 
the energy laws and discusses what quantities are preserved with the numerical schemes 
being used. Numerical results for certain bubble test problems and an MHD simulation 
are presented in section [5j where it is shown that this approach is accomplishing the 
intended goal. Finally, a discussion of the results and summary is given in section [6] 

2. The Two-Phase Equations 

The dynamics of mixing between two fluids plays an important role in many physical 
applications. Here, we use a model that couples the Navier-Stokes equations with a phase- 
field transport equation. We use the Allen-Cahn-type equations, which were originally 
introduced in [13] to describe the motion of anti-phase boundaries in crystalline solids. 
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The Allen-Cahn equations have since been used in many interface problems in fluid 
dynamics through a phase- field approach [31 SI El 031 US] ■ They are 



^ + (u-V)u + Vp + AV-(V0<8)V0)- ) uV 2 u = 0, (f) 

V-u = 0, (2) 
^+u.V0- 7 (v 2 ^i0(0 2 -l)) = 0. (3) 



Here, u is the fluid velocity, p is the fluid pressure, and <j> is the phase-field function. The 
physical parameters are the mixing energy parameter, A, the viscosity of the two fluids, 
/i, the elastic relaxation time of the system, 7, and the mixing layer width, e. We assume 
that both fluids have the same viscosity and mass density (which is assumed to be 1). It 
should be noted that as e — > 0, the parameter A is associated with the surface tension at 
the boundary. 

With the "mixing" energy, W(<f>), described above, and following the general approach 
of Onsager (TH |T71 [T51 US] , we can prescribe the whole system as a dissipative system: 

I H l -p{<P)\u\ 2 + \Wm^ = -^ M (0)|Vu| 2 + ^|0| 2 dx, 

where |0| 2 dx represents the microscopic internal damping, with the kinematic as- 

d(j) 

sumption of <fi — — — h u • V0. We can then employ the least action principle (LAP) 
at 

for the Hamiltonian part of the system and the maximum dissipation principle (MDP) 
for the dissipative part of the system (T6l [T71 [TBI 03 EQ] . In the case that \i\ = [12 = /i 
and pi = 1, we obtain equations (jl}-(|3|. In case 7 = 0, i.e., we neglect the internal 
damping, we can show that ([l])-([3]) approaches the sharp interface formulation [T51 |2T] . 
In particular, 

AV • (V0 ® V</>) = A (\7 2 (f>V<j> + V 

where A (V 2 </>V0) reflects the surface tension, fcf/n, modulating the pressure contribu- 
tions. 

Remark Later, we indicate the structural similarities between the diffusive interface for- 
mulation and that of resistive MHD systems. In this setting, the terms AV • (V<^> ® V</>) — 

V J are equivalent to the Maxwell stress, and AV 2 0V0 is equivalent to the Lorentz 

force. 

Next, we discuss the numerical methods that are used to solve this problem. 



3. FOSLS Formulation 

To solve the nonlinear system of equations, which for discussion is denoted by C(u) = 
/, a first-order system least-squares discretization (FOSLS) is used [5] [52]. First, 
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consider a linear first-order system, denoted by Lu — /. Then, the linear problem is 
recast as the minimization of a functional constructed by taking the L 2 norm of the 
residual of each equation. This is written as 

u* = argminG(u; /) := argmin \\Lu - (4) 
where is the solution in an appropriate Hilbert space V. 

Thus, u* satisfies [G'(u*)]v = 0, which is the Frechet derivative of G at applied to 
V G V. This results in the weak form of the problem: 

find Uit € V such that 

(Lu*,Lv) = (f,Lv) VweV, (5) 

where (•,•) is the usual L 2 inner product. Desirable properties of the bilinear form, 
(Lu, Lv), are 

continuity (Lu,Lv) < C2||u||y|Hlv Vn,«e V, (6) 

coercivity (Lu,Lu) > ci||u||y Vu£ V. (7) 



These properties imply the existence of a unique solution, », £V, for the weak problem 
©■ 

Next, we approximate it* by restricting (Wj) to a finite-dimensional subspace, V h C V, 
which leads to (5.) restricted to V'\ Since V" is a subspace of V, the discrete problem 
is also well-posed. Choosing an appropriate basis, V h — span{$j}, yields an algebraic 
system of equations involving the matrix A with elements 

(A)ij = (L* 3 -,i*j). 

In general, if V is a subspace of H 1 (product space), and continuity and coercivity hold 
we say that the FOSLS functional is H 1 equivalent. In this case, the linear system is 
amenable to an iterative solution by multilevel techniques. In particular, AMG has been 
shown to work well on a wide range of such problems [51 171 151 [51 125 1 135 1 1^1 12711251 125] . 

In addition, note that the functional yields an a posteriori error measure. If u h G V' 1 , 
then 

G(u h ;f) = \\Lu h -f\\ 2 = \\Lu h - Lu.A\ 2 = \\Le h \\ 2 a c\\e h \\ 2 v . (8) 

The last relation comes from the continuity and coercivity bounds in |6]) and (|7|). Thus, 
the functional value is equivalent to the error measured in the Hilbert space norm. In 
general, the constant c in equation ^ depends on the continuity and coercivity constants, 
C2 and c\. These constants, of course, depend on properties of the PDE as well as on 
boundary conditions and the computational domain. For this paper, we study simple 
domains and boundary conditions for which these constants are not large. In addition, 
our numerical results show that the problem parameters do not adversely affect our 
accuracy either. Addressing these issues analytically would be a topic of another paper 
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and has been studied in other problems such as in [51 H3 HO] • Locally, on any subset 
of fi, such as an element, the functional yields an estimate of the error. This property 
of FOSLS helps make it possible to efficiently solve complex systems. At each step in 
the solution algorithm, a local measure of the functional is available. Since this is the 
norm we are minimizing, this allows judgements to be made based on estimates of the 
increase of accuracy that results from an increase in computational cost. As a result, 
methods such as nested iteration and adaptive local refinement can be used to improve 
the efficiency of solving such systems. This has been applied to various complex fluids 
problems, including magnetohydrodynamics [TTJ [T^l [301 [3T] . A brief description of the 
algorithms used is described in section [5] 

3.1. Linearization (Newton-FOSLS) 

Since the two-phase system is nonlinear, we first linearize it before we put it into a 
FOSLS weak form. The system becomes 

£(ti +u)«£(u ) + [£'(V)R (9) 

where Uq is the current approximation and [£'(uq)] is the Frechet derivative evaluated 
at uq. The functional then becomes 

G(uJ) = \\[C'(u Q )]u-(f-C(u Q ))\\ 2 . (10) 

Minimization of the linearized functional yields that satisfies the weak form: 

find i, 6 V such that, 

([C'(u )}u t ,[C'(u )}v) = {(f-L(u )),[C'(u )}v) VveV. (11) 

Once u is found, it is added to the previous iterate to get the next approximation, 

til = u o + u. (12) 

Thus, in a nonlinear setting, the FOSLS approach can be applied and, if the linearized 
functional retains continuity and coercivity, ([6| and ([7]), it also retains the desirable 
properties as described for linear systems. In addition, the nonlinear functional, that is, 
the L 2 norm of the residual of the nonlinear system, can be computed as well as the linear 
functional or the L 2 norm of the linearized system. This allows the two functional values 
to be compared after each linearization and helps determine if the Newton iterations are 
converging as expected. 

A similar approach would be to create the FOSLS functional of the nonlinear problem 
and linearize this functional instead. This FOSLS-Newton approach generally involves 
more terms than the Newton-FOSLS method described here. As the approximation 
approaches the solution, these additional terms are higher-order and tend to zero faster 
than the overall functional and the two methods tend to be the same. Thus, when 
nested iteration is used, as it is here, there tends to be very little difference in the 
convergence behavior of these two approaches. FOSLS-Newton may be more robust in 
some especially challenging situations, but the Newton-FOSLS approach is simpler and 
has been successful in a number of applications [251 32 , so we confine our presentation 
to it in this paper. 

6 



3.2. Velocity- Grad Formulation 

To reformulate 0-([3]) as a first-order system, we introduce new variables that repre- 
sent the gradient of the velocity vector, u, and the phase-field function, <fi: 



V = Vu 



"21 V 2 2 

B = V</> = 



d x ui d x u 2 

dyUi 8yU2 



d x 4> 
d.,6 



(13) 



(14) 



This formulation has been shown to work well for the Navier-Stokes equations alone 
when velocity boundary conditions are given 8,9. To make the system H 1 equivalent, 
a few auxiliary equations are added to the system. These are consistent with the original 
equations and make the first-order system well-posed. The resulting system is as follows: 



at 



du 

~dt 



u B 



V- Vu 


= 2> 


(15) 


V x V 


= 0, 


(16) 


V ■ u 


- 0, 


(17) 


V (trV) 


= 0, 


(18) 


Vp + AB (V • B) - fiV ■ V 


= 0, 


(19) 


B — V0 


= 0, 


(20) 


V x B 


= 0, 


(21) 


7(v-B-i^(0 2 -l)) 


= 0. 


(22) 



In practice, equation ( 18 ) is achieved by eliminating one of the diagonal components 
of V_ in terms of the other diagonal components. This comes from the divergence- free 
constraint on the velocity. Thus, in 2D, there are 13 equations and 9 unknowns and, in 
3D, there are 29 equations and 16 unknowns. 

3. 3. Equivalence of the FOSLS Functional 

As described above, if continuity, ([6]), and coercivity, Q, of the FOSLS functional 
holds, we say the functional is H 1 equivalent. This can be shown for the linearized system 
at each time step. We assume that an implicit backward difference time-stepping scheme 
is used and consider this semi-discrete system first. This says that, for each Newton step 
and time step, the FOSLS discretization gives a well-posed numerical system to solve. 
We state the main theorem here and give the proof. However, many of the details are 
omitted to save space and since they have been shown in previous works. 

Theorem 3.1 (cf. [5J [7J |H1 M E3] ) • Let C be the linearized first order two-phase flow 
operator as in equations fifty - ffBOfy and let G be the FOSLS functional as defined in equa- 
tion (10) of this linearized system. Let the domain, il, be a bounded convex polyhedron 
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with connected boundary or a simply connected region with a C ' boundary. Denote 
U = (u, y, p, (j), B) T and let 

uev-.= [h 1 ^)] 2 x [H 1 ^)] 4 x [ig(n) n fr x (n)] x [ir^(n)] x [h 1 ^)]. 

Assume u™, F™, B™, and 0" are given members of our finite- element space and, there- 
fore, 

ll u ™llooj ||Y_ n ||oo> | IB™ | |oo , ||V • B^loo, and \\<p n \\oo < oo. Then there exists positive 
constants, C\ and c%, that yield the following bounds: 



a. 



b. 



G< C2 (||u||? + ||Z||? 



|B||? 



I), 



(23) 



ci(||u||f + ||Zllf + IWIf + l|B||f + ||0||f)<G. (24) 

These bounds mean that the functional, G, is equivalent to the H 1 norm of the error. 
The spaces above are defined as follows: 



Lg(fi) = {pe L 2 (n) : [ pdx = 0}, 



H^Q) = {pei 2 (0):V P e(L 2 (0)) 3 }, 
H 1 (f2) rf = the d-dimensional version of i? 1 (f2) , 
Ho(n) = {p G L 2 (n) : Vp G (L 2 (n)f and p = on dft}. 



Proof of\S.l\ Proving continuity, ( 23 ), of the linearized first-order system involves several 



instances of the triangle inequality and is, therefore, easy to show. 



Coercivity, (24 1, is accomplished by treating the two coupled block systems, Navier- 



Stokes and the phase field equation, separately. A linearized system of (15)-(22) can be 
written in block form as follows: 



MA 



Ans Anp 
Apn A p h ase 



I u \ 

p 



V B J 

Here, Ans represents the linearized Navier-Stokes equations in first-order form and with- 
out the coupling term, ( 15 )- ( 19 1 . Similarly, A p h a se represents the linearized phase func- 



tion equations without the velocity field coupling, (20 1- (22). The coupling terms are as 
follows: 



A NP B = AB"V • B + A(V • B")B, 
A PN u = B" u. 
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The superscript n represents that the solution is from a previous Newton step. Since we 
are working on discrete finite element spaces, we can assume the L°° norms of all the 
previous iterations are bounded. 

In 0H], coercivity has been shown for the A^s block given here. For the phase field 
system, the coercivity bound can be shown in a similar fashion. In fact, ignoring the 
nonlinear terms, the system resembles the FOSLS formulation used on the heat equation 

Each 

B — = 0, 
V x B = 0, 
= on c?f2, 
n x B = on dfl. 

Forming the functional of the linearized system for just the phase field block, we write 

G ph ase = |K + /3 • B - 7 V • B||g + ||V x B||g + ||B - V0||g. 

Again, it is assumed that an implicit time-stepping scheme is used so a depends on the 
time step. With part of the nonlinear coupling included in this block, (3 depends on the 
parameters of the problem and the previous solutions. This in fact makes the problem 
resemble an advection-diffusion type system. Since this functional only adds zeroth-order 
terms (i.e. no derivatives) to this system, a standard compactness argument, [34], can 
be used to show that there exists a constant, C > 0, such that 

C(||B||? + ||0||?) <G phase . 

It should be noted that this requires c||B||f < ||V ■ B||q + ||V x B||q, for some constant 
c > 0. However, this is true due to the assumptions on the boundary and the boundary 
condition, n x B = [33]. Thus, the uncoupled system is H 1 equivalent: there exists 
positive constants, cl and cjj, such that 

2^-ll/ A-NS \ 7ju2 » \\i,\\Vt 



c L \\u\\i<\\^»* A ^ ase jUW^cuWUWf 

Next, we consider adding the off-diagonal coupling blocks of the system. In [53], the 
following lemma was proved: 

Lemma 3.2. Let C be a 2 x 2 upper triangular block matrix such that 

uch: 

Let Ai and A 2 be invertible and let there exist a positive constant, C , such that 

\\Tn 2 \\l<C\\A 2 n 2 \\l 



Then 





1 




o 1 


A 2 J 





l2 <(l + C)\\< Al T ^ ' " ; 



A 2 J \ u 2 



o- 



Looking at the upper triangular portion of the two-phase flow system, we see that 



\\A NP B\\ 




AB"V ■ B + A(V ■ B n )B\\ 2 Q < Ci||B||§ + C 2 ||V • B 



for some positive constants C\ and Ci. Adding a curl term at the end and using a 
Poincare inequality, it is easily shown from the coercivity of the phase-field block that 



Finally, to show equivalence of the full linearized system, it should be noted that 
the lower left block, Apjy, involves a term with no derivatives and another standard 
compactness argument can be used. Here, again, we must assume that the previous 
iterations of u and B are bounded. Thus, there exists positive constants c\ and c 2 such 
that 



Therefore, for each time step and linearization, there exists a weak solution to our 
linear system, and our multilevel solvers will converge to that solution. The constants 
ci and C2 depend on physical parameters such as 7, fi, and e, as well as the previous 
Newton step iterations and time step approximations. However, they do not depend on 
the number of degrees of freedom or the size of the finite element grids being used. In 
the next section, we discuss the energetics of complex fluid flow. 

4. Energetics 

While most physical systems can be described by a set of PDEs, usually an underlying 
energy law also describes the system. In many cases, the PDE is actually derived from 
the energy laws via a calculus of variation (see e.g. [2Q])- The main idea is that the 
change in the total energy of a system over time must equal the total dissipation of 
the system. If the system is conservative, or there is no dissipation, then the total 
change is zero and it is a Hamiltonian system. Thus, the interaction or coupling between 
different scales or phases plays a crucial role in understanding complex fluids. Any set 
of equations that describe the system should then be a result of the energy laws. The 
energetic variational approach (EVA) takes the energy laws of a complex fluid model and, 
using the least action principle (LAP) and the maximum/minimum dissipation principle 
(MDP) [Tni Q~71 [TS1 [05], yields a weak form of the system, which we approximately 
preserve in our numerical simulations. The energy variation is based on the following 
energy dissipation law for the whole coupled system: 






□ 



Qj^total 



V, 



(25) 



dt 
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where E total is the total energy of the system and V is the dissipation. The LAP, which 
is also referred to as the Hamiltonian principle, or principle of virtual work, gives us 
the Hamiltonian (reversible) part of the system related to the conservative force. At the 
same time, the MDP gives the dissipative (irreversible) part of the system related to the 
dissipative force. 

In [H [3 [13j HU [15] and others, an energy law is developed for the Allen-Cahn phase- 
field model described above. The total energy is a combination of the kinetic energy 
(driven by the fluid) and the internal energy (driven by the interface): 



E 



total 



l 



U" 



XW{4>) \ dx, 



where W(<p) is the Ginzburg Landau mixing energy described above [1], thus representing 
a competition between two fluids with their hydrophilic and hydrophobic properties. The 
dissipation is a result of a diffusive term from the viscosity of the fluid and from the 
diffusion of the interface itself: 



V = 



M |Vu| 2 + -|0| 2 ldx = 

7 J Jn 



//|Vu| 2 + A 7 



V 2 4>- 



1 



-1 



dx. 



Thus, any system describing this behavior, such as equations 0-([3| or the FOSLS 
formulation ( 15 )-( 22 ) , should globally produce results that approximately satisfy the 



following energy law: 



d 
dt 



u| 2 + XW{cf>) Ux = - / { ^|Vu| 2 + A 7 



(26) 



The overall energy law reflects the multi-scale, multi-physics nature of the system. As 
mentioned above, the variable <$> represents the microscopic description of the mixture. 

dej) 

The kinematic assumption, 6 — — — H u • V6, stands for the influence of macroscopic 

dt 

dynamics on the microscopic scale. The form of the internal dissipation is the consequence 
that we are looking at the long (or macroscopic) time dynamics. This brings us to the 
equation for tfi, 



90 
di 



u • V(f> = 7 V 



1 



1 



This "gradient flow" dynamics is another formulation of the near equilibrium, linear 
response theory [T5J [TH] . 

As for the macroscopic force balance, the LAP gives the following Hamiltonian system, 
with 7 = 0: 



i:)u 

at 



(u-V)u + Vp = -AV • (V0 ® V<£) , 

V • u = 0, 

d 4 + n-V4> = 0. 
at 
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This describes the transient dynamics of the system. 

On the other hand, the MDP yields the dissipation of system, which stands for the 
long macroscopic time dynamics: 



-/iV 2 u + Vp 



A . 

7 



at 



v ■ u = o, 

u-V0 = 7 (V 2 0+4(0 2 -l) 



Remark ±<j>V<P = A (V 2 + £ (</> 2 - l) 0) V</> = AV 2 0V</> + AV (± {(f) 2 - l) 2 ) . This 
shows the consistency between the LAP and the MDP 0]. 

System ([l])-([3| is really the hybrid of these two systems. 

Remark It is worth pointing out that due to the generality of the formulation, analo- 
gies of (Tj)-([3) can be found in many different physics. In particular, without the 
In ie 1 — l) 2 ^ x term, the system is equivalent to the resistive MHD system. In fact, 
it is this understanding that allows us to arrive at a similar weak form for the FOSLS 
algorithms as in [TT]. The choice of using the gradient of <j> as an auxiliary variable in 



( 14 ) was made for this reason. 



Finally, we point out that the advantage of the above energetic variation approach 
enables us to derive the thermodynamic-consistent models involving different physics at 
different scales. Restricting to a finite-dimensional space, we show numerically that the 
FOSLS discretization method is capable of globally preserving this energy law. 



5. Numerical results 

In [THE], an algorithm is devised to solve a system of nonlinear equations, C(u) = /. 
Starting on a coarse grid with an initial guess there, the system is linearized and the 
linearized FOSLS functional is then minimized on a finite clement space. At this point, 
several algebraic multigrid (AMG) V-cycles are performed until there is little to gain in 
accuracy-per-computational-cost. The system is then relinearized and the minimum of 
the new linearized FOSLS functional is searched for in the same manner. After each set 
of linear solves, the relative difference between the computed linearized functional and 
the functional of the nonlinear system is checked. If they are close and near the minimum 
of the linearized functional, then it is concluded that we are close enough to the minimum 
of the functional of the nonlinear operator and, hence, we have a good approximation to 
the solution on the given grid. Next, the approximation is interpolated to a finer grid 
and the problem is solved on that grid. At this stage, the grid can be refined uniformly 
or locally. This process is repeated to yet finer grids until an acceptable error has been 
reached, or until we have run out of computational resources, such as memory. If, as 
in the case of the 2-phase flow equations, it is a time-dependent problem, the whole 
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process is performed at each time step. This algorithm is summarized in the following 
pseudocode. 

Algorithm 1: Nested Iteration Newton-FOSLS AMG 

for t = 1 to max time step do 
Go to coarsest grid 

while fine grid resolution is not reached do 
for n — 1 to max Newton step do 
Linearize first-order system 
Minimize FOSLS functional 

while functional of current approximation is large do 
I Solve Ax = b with AMG 
end 

Compute relative difference between linearized and nonlinear operator 
if small then 
| Exit Newton loop 
end 
end 

Refine grid 
end 

Update time step 
end 

The main goal behind this algorithm is to reduce the amount of work needed to solve 
the system of equations. Relatively cheap work is done on the coarse computational grid 
to get better starting guesses for the solution of the system on the finer finite element 
grids. Then, less iterations are needed at the resolution that would require the most 
work. The algorithm described above also takes into account the fact that resolving too 
much on coarse grids is unnecessary. Using the FOSLS functionals as estimates of the 
errors, stopping criteria is implemented for each stage of the algorithm, optimizing the 
amount of work per computational cost. 

In addition to the nested iteration algorithm, adaptive local refinement is used. Again, 
using the FOSLS a posteriori error estimates, an efficiency-based adaptive local refine- 
ment scheme, known as ACE (30J, [35] can be used fairly easily. The idea behind this 
scheme is to predict the amount of work and the error reduction obtained from refining 
only certain portions of the computational domain. Then, the optimal grid is chosen. In 
simulations of magnetohydrodynamics, which have similar structure to the complex fluid 
problems described here, this scheme has been shown to reduce the amount of work by 
a factor of ten in the FOSLS framework [3TH [3Tj . 

For the following test problems, we use the FOSLS finite element method along with 
nested iteration, AMG, and adaptive local refinement. Here we show that, using these 
methods, we are able to solve the complex fluid problems well, while still adequately 
preserving the energetics of the system. 

5.1. Two-Phase Flow 
5.1.1. Coalescence 

In this test problem, we start with two bubbles of a certain phase immersed in a fluid 
of a different phase on the unit square, f2 = [0, 1] x [0, 1]. The initial interface between 
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the two fluids is given such that the two bubbles are osculating or "kissing", and then 
the system is allowed to evolve. As a result, the two bubbles begin to coalesce into one 
and, eventually, due to the dissipation in the system, the newly formed bubble slowly 
shrinks and is absorbed by the outside fluid. The initial conditions for this system are 
similar to that in [4] and are as follows: 



= tanh 



di(x,y) 
2 V 



tanh 



2 V 



-1, 



where 

di(x,y) = y/{x - 0.38) 2 + (y - 0.5) 2 - 0.11, 

and 

d 2 (x,y) = y/(x - 0.62) 2 + (y - 0.5) 2 - 0.11. 

The parameters used are rj = 0.01, fi = l,e = 0.01, 7 = 0.01, and A = 0.0001. 

A second-order fully implicit backward differencing formula (BDF-2) is used and 
allows us to take a larger time step for the simulation than in previous work [4] . For the 
test problem shown, 100 time steps of size At = 0.01 are taken. Plots of the phase field, 
0, are given in figure [2] Here, one can see the coalescence of the two bubbles over time 
and the eventual dissipation. In addition, more refinement is done in the region of the 
interface. Therefore, most of the computation is being focussed on resolving this region, 
which is driving the physics in the system. 

t = 0.01 t = 0.1 t = 0.2 



OO CO CO 



t = 0.4 t = 0.7 t = 1.0 



CD O o 



Figure 2: Phase field, <j>, at time t = 0.01, 0.1, 0.2, 0.4, 0.7, and 1.0. Red represents 1.0 and blue 
represents -1.0. 

To quantify how well the nested iteration algorithm with adaptive local refinement 
performs, the simulation was run using uniform refinement. An estimate of the work is 
computed by calculating the number of non-zeroes in each discrete system and determin- 
ing how many equivalent fine-grid matrix-vector operations are performed. We call this 
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measure a work unit (WU), and it is a measure of the cost of performing one relaxation 
sweep, such as Gauss-Seidel or Jacobi, on the finest grid of the simulation. Table[T]shows 
that, while using adaptive refinement requires more work units (i.e. more V-cycles are 
needed to solve the problem), on average the number of non-zeroes on the finest grid 
is about one-tenth of that for uniform refinement. This corresponds to a work ratio of 
about 0.21. In other words, using adaptive local refinement saved about 80% of the com- 
putational cost. Therefore, on average, the above strategies on this test problem cost the 
equivalent of performing less than 49 Gauss-Seidel relaxation sweeps on a quadratic uni- 
form grid with 65, 536 (256 x 256) elements. It should be noted here that the functional or 
estimate of the error for adaptive refinement is about half of the functional obtained from 
using uniform refinement on average. Both schemes got below the prescribed tolerance, 
but, on many time steps, the adaptive algorithm got to a much lower error (occasionally 
by an order of magnitude). Therefore, as far as accuracy-per-computational cost, the 
adaptive refinement scheme may be even better than the numbers indicate. What we 
can say is that it used less than 9% of the elements and took almost one-fifth of the 
computational time to solve the problem. 



Avg WU 


Avg Nonzeros 


Avg Functional 




Uniform 




227.1 


281,292,867 


0.456 




Adaptive 




530.56 


24,955,707 


0.265 


Work Ratio 
0.213 




Element Ratio 
0.089 



Table 1: Average estimates of work and accuracy using uniform refinement versus adaptive local refine- 
ment 

In addition, Table [2] shows that, for each refinement level, the number of Newton 
steps needed to solve the nonlinear problem decreases from almost 5 on the coarse grids 
to 1 on the finer grids. For all time steps, on the finest-grid, only one Newton step was 
needed. This greatly reduces the cost of the algorithm compared to if one were to only 
solve on the finest grid. 

Next, we see in figure [3] that the energy law of this system is satisfied during the 
simulation. The total conservative energy over the whole domain is decreasing and the 
change in energy coincides with the dissipation in the system. Thus, the numerical 
simulation is resolving the expected physics of the system. 

5.1.2. Square Bubble 

In this test problem, we start with a single bubble immersed in a fluid of a different 
phase on the unit square, ft = [0, 1] x [0, 1]. The initial interface between the two fluids 
is given such that the bubble is in the shape of a square. As time evolves, the interface 
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Grid 


Elements 


Avg Newton Steps 


10 


4 


4.6 


9 


4 


2.2 


8 


16 


2.9 


7 


47.1 


2.7 


6 


84.9 


2.3 


5 


178.1 


2.0 


4 


417.7 


2.0 


3 


1,005.0 


1.1 


2 


2,084.9 


1.0 


1 


5,812.4 


1.0 



Table 2: Average Newton steps per grid level over all time steps. Level = 1 is the finest grid. The first 
grid (level 10) uses bilinear elements. 




Time 

Figure 3: Energy laws for coalescence problem. Top: Total conservative energy over time. Bottom: 
Change in energy (asterix) over time versus negative dissipation (circle) 



forces cause the bubble to transform to a circular shape before diffusing away. In the 
process, oscillations between a square and diamond shape are formed. This reflects the 
elastic short-time behavior of the system. The initial conditions for this system are 
similar to that in [2] and are as follows: 





/ 

+1 (x, y )e[-i,i]x[-i,i] 

— 1 otherwise 



The parameters used are fi = 0.1, e = 0.02, 7 = 0.01, and A = 0.1. 

Again, a second-order fully implicit backward differencing formula (BDF-2) is used. 
For the test problem shown, 100 time steps of size At = 0.01 are taken. The phase 
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field is shown in figure [4] As in the previous example, the refinement is being placed in 
the appropriate regions and requires about 20% of the work that would be required if 
uniform refinement were used. 



t - 0.01 t = 0.03 t = 0.07 




t - 0.1 t = 0.18 t = 0.25 



o o o 



Figure 4: Phase field, cj>, at time t = 0.01, 0.03, 0.07, 0.1, 0.18, and 0.25. Red represents 1.0 and blue 
represents -1.0. 

Due to the bigger jump at the interface, though, more Newton steps are needed in 
this problem than in the first. However, using nested iteration results in many of these 
linearizations being performed on the coarser grids. As an example, at time step 3, when 
the bubble begins to lose its square shape, the number of Newton steps goes from 10 to 
1 as we move up through the grids. See table [3] 



Grid 


Elements 


Newton Steps 


9 


16 


10 


8 


64 


4 


7 


112 


5 


6 


208 


5 


5 


463 


5 


4 


1,423 


4 


3 


3,766 


3 


2 


12,793 


3 


1 


40,030 


1 



Table 3: Newton steps per grid level for time step 3. Level = 1 is the finest grid. 

In figure [5] one can see that the energy laws are being satisfied once again. An initial 
increase in energy occurs before the dissipation in the system takes over, causing the 
energy to decay. This is most likely due to the fact that the initial condition of the phase 
field is discontinuous at the interface. Once the simulation is started, this function is 
immediately smoothed out and the energy behaves as expected. By increasing 7, the 
elastic relaxation time in the system, the dissipation is increased and actually causes the 
bubble to dissipate away much faster. In fact, the entire bubble vanishes after time step 
30 and the total energy goes to zero. This can be seen in the energy plots in figure [6] 
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The tiny jump in energy appears to occur when the region containing the bubble finally 
disappears and the system is changed to that of one fluid. With a smaller time step, this 
transition occurs more smoothly. 



Total Conservative Energy: ia = 0.1 y=0.01 \ = 0A, BDF2dt = 0.01 



800 
600 
400 
200 


-200 
-400 



Change in Energy over Time and Dissipation 



-dE/dt 
--Disp 



Figure 5: Energy laws for square problem (7 = 0.01). Top: Total conservative energy over time. Bottom: 
Change in energy (asterix) over time versus negative dissipation (circle) 



Total Conservative Energy: \i = 0.1 y = 0.1 % = 0.1 , BDF2 dt = 0.01 



500 


-500 
-1000 
-1500 



0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 

Time 

Change in Energy over Time and Dissipation 




0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 

Time 



Figure 6: Energy laws for square problem (7 = 0.1). Top: Total conservative energy over time. Bottom: 
Change in energy (asterix) over time versus negative dissipation (circle) 



5.2. MHD 

As described above, the characteristics of the two-phase system is similar to that 
of resistive MHD. The main results and description of this problem using the FOSLS 
and nested iteration framework can be found in [TTJ [T2J |3D] . However, here, we briefly 
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mention the energetics of the system. For this system, the energy law must satisfy the 
following equation: 

^/ {iM- + i|B|'}*-/ o {^> + ^|VB|.} fc 

Here, u represents the fluid velocity and B represents the magnetic field. R e and Sl are 
the fluid Reynolds number and magnetic Lundquist number, respectively. 

A magnetic reconnection problem is simulated using the above numerical algorithm. 
Figure [7] shows that the discrete energy law is being satisfied. The total conservative 
energy decays as a direct result of the dissipation in the system. The oscillations in the 
change in energy plot are due to a lower order differentiation being used to compute 
after the simulation was complete. Similar to the first two examples, it took on average 
of 10 WU to solve the problem using nested iteration and adaptive local refinement. 



Total Conservative Energy: R = 50001 S = 50001 , BDF2 dt = 0.1 
1.696 1 1 1 1 1 1 1 




01 23456789 10 

Time 

v in" 3 Change in Energy over Time and Dissipation 




Figure 7: Energy laws for mhd test problem. Top: Total conservative energy over time. Bottom: Change 
in energy (asterix) over time versus negative dissipation (circle) 



6. Discussion 

This paper shows that using the Newton-FOSLS finite-element formulation, along 
with nested iteration, algebraic multigrid, and adaptive local refinement, accurately re- 
solves two-phase flow simulations while still preserving the energetics of the system. The 
numerical algorithms used were designed in a general setting and have the advantage 
of reducing the amount of computational cost. For the two-phase flow problems, it is 
possible to resolve the physics in the equivalent of less than 50 fine-grid relaxation sweeps 
for each time step. This is accomplished by reducing the amount of work done on highly 
resolved grids and by reducing the number of elements needed to solve the problem via 
adaptive local refinement. As a result, complicated systems of PDEs can be solved very 
efficiently without losing the ability to measure certain physical quantities. This paper 
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shows that, when applied to two-phase flow simulations, all the advantages of using im- 
plicit finite element methods are obtained, and the FOSLS methodology is capable of 
capturing the important aspects of the model, such as the energy laws. 

Despite the success of the above algorithms in capturing the correct energy laws effi- 
ciently, there has been one bottleneck. In all the simulations, classical algebraic multigrid 
(AMG) was used to solve the discrete linear systems. Classical AMG was designed for 
scalar elliptic problems and M-matrices. As can be seen in figure |5J the average conver- 
gence factors for AMG over each time step can be quite poor when it is used on systems 
for which it was not designed. While it is always bounded away from 1, one can see that 
the worse the AMG convergence factor gets, the more computational work is needed to 
solve the problem. Therefore, future work will involve investigating various multigrid 
methods that will be more well-suited to the test problems described above. Algebraic 
methods that are designed more specifically for systems of PDEs will be studied as well 
as geometric multigrid. All the domains have been well structured and, therefore, it 
is possible to use more of the geometric information to create a better solver. If it is 
possible to improve the linear solver using these methods, then the amount of work can 
be reduced even further. 




0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 



Figure 8: Work units per time step as well as average convergence factor per time step. 
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